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Abstract. We suggest a novel shape matching 
algorithm for three-dimensional surface meshes of 
disk or sphere topology. The method is based on 
the physical theory of nonlinear elasticity and can 
hence handle large rotations and deformations. De¬ 
formation boundary conditions that supplement the 
underlying equations are usually unknown. Given 
an initial guess, these are optimized such that the 
mechanical boundary forces that are responsible for 
the deformation are of a simple nature. We show 
a heuristic way to approximate the nonlinear opti¬ 
mization problem by a sequence of convex problems 
using finite elements. The deformation cost, i.e., the 
forces, is measured on a coarse scale while ICP-like 
matching is done on the fine scale. We demonstrate 
the plausibility of our algorithm on examples taken 
from different datasets. 

1 Introduction 

1.1 Motivation 

Shape matching is a difficult but nonetheless 
important problem in computer vision, medical 


imaging, and computer graphics. It is the key 
ingredient for recognition, retrieval, alignment 
of scanned data, information transfer, shape in¬ 
terpolation, statistical shape modeling, space- 
time reconstruction and more. 

This work aims to contribute a novel frame¬ 
work to solve for unknown shape deformations 
of a pair of two-dimensional surfaces (source and 
target) in three dimensions. Our method is built 
on the observation that in applications surfaces 
often represent the boundaries of actual physi¬ 
cal entities, i.e., surfaces of elastic bodies. Shape 
change can therefore be explained by means of 
forces acting on the specific elastic body. Their 
magnitude can be interpreted as a measure of 
how “difficult” it is to achieve a certain shape 
change. Elastic models are, we believe, well 
suited to give valuable information about shape 
changes among different objects. In particular, 
we believe, that sparsity of forces is a good prior 
for explaining an observed deformation within a 
semantic class of objects, e.g., if we were to com¬ 
pare two different shape articulations. 

From a mathematical point of view shape 
matching is an ill-posed inverse problem and is 
usually tackled in a heuristic manner. It is com¬ 
putationally challenging since the search space 
for correspondences is usually large. Hence, 
there is a need to explore this search space in 
a reasonable manner. 

1.2 Our Method 

In |3S] the authors intuitively describe the de- 
sired shape blending algorithm as one which re¬ 
quires the least work to deform through bend¬ 
ing and stretching. Many shape matching algo¬ 
rithms are designed in this spirit as one usually 
seeks “maximum alignment by means of mini¬ 
mal cost” in an optimization framework. Our 
work uses the physical theory of nonlinear elas¬ 
ticity and is based on the same principle as will 
be elaborated in the following. 

In this work we suggest a method for align¬ 
ing two surface meshes (manifolds) embedded in 
three-dimensional space. We use the observa¬ 
tion that surface meshes often represent (parts 
of) the boundaries of actual physical bodies. 
Hence, the cost of a deformation of a surface 
can be measured by forces acting on the bound¬ 
ary of the solid body described by the surface. 

We believe that in many realistic scenarios 
even complicated deformations take place due 
to or can be explained by means of surpris- 
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ingly simple forces. Examples are the change 
of a pose of an animal or a human where forces 
act mainly only on the articulated parts, i.e., 
the acting forces can be regarded as sparse. We 
furthermore anticipate that, given two semanti¬ 
cally similar objects, seeking a deformation that 
is subject to sparse and isotropic surface forces 
is a reasonable strategy to match them. These 
forces could then be used to intuitively compare 
two different deformations and to assess their 
severity. 

The underlying equations of nonlinear elas¬ 
ticity that we employ form a system of nonlin¬ 
ear partial differential equations (PDEs) and are 
subject to boundary conditions. These can be 
given as deformations and/or forces prescribed 
on the boundary of the body. Given two shapes 
boundary conditions of either type are usually 
unknown but one can often give an initial guess 
of boundary correspondences, i.e., of the bound¬ 
ary deformation. Such an initial guess can then 
be optimized according to our prior on the forces 
to obtain a “better” set of boundary correspon¬ 
dences that supplement the PDE. The initial 
guess plays the role of a spring force between 
source and target shape and a penalty needs to 
be paid for deviation. The optimized bound¬ 
ary correspondences can then be used to solve 
the underlying elasticity PDE using a common 
Newton-Raphson scheme to update the full vol¬ 
umetric deformation of the body and the entire 
procedure can be iterated until a good match is 
found. 

Using the popular finite element method 
(EEM) we show how to connect the boundary 
deformation to the boundary forces by means 
of a splitting of the tangent stiffness matrix, lo¬ 
calized around a certain deformation. This so- 
called condensation is necessary to reduce the 
optimization to the boundary of the body and 
has the advantage that it lowers the dimen¬ 
sion of the optimization significantly. Eurther- 
more, the condensation allows us to formulate 
the problem as a sequence of (convex) second 
order cone problems. 

As elaborated above we aim to match surfaces 
represented by triangular meshes but the EEM 
works on volumetric meshes. Hence, we create 
a tetrahedral mesh of the body represented by 
the triangular mesh. The surface meshes that 
we match in our examples have between 9000 
to 55000 triangles and constructing a tetrahe¬ 
dral mesh that resembles the surface density of 
meshes would result in a computationally very 


large problem. Therefore, we create a coarse 
tetrahedral approximation of the surface. Every 
point on the surface mesh can then be coupled 
to the coarse boundary mesh of the volumetric 
mesh. A deformation of the coarse boundary 
then induces a deformation of the surface. This 
way we can prescribe a matching (in particu¬ 
lar the initial guess) on the surface mesh while 
measuring the deformation cost on the coarser 
volumetric mesh. 

1.3 Previous Work 

1.3.1 Overview 

An optimization based shape matching algo¬ 
rithm usually consists of two ingredients. Eirst, 
one needs a regularizer that determines (the 
preference of) the class of deformations in con¬ 
sideration. Well investigated is the field of rigid 
matching, i.e., the admissible transformations 
are translations and rotations [Il|25l|39]. This 
is a common scenario if one is given a range 
scan of rigid objects that can be aligned with 
a single rigid transform. Eor non-rigid match¬ 
ing common regularizers that are used in the 
literature are the popular As-Rigid-As-Possible 
(ARAP) functional [2] and the As-Afhne-As- 
Possible (AAAP) functional [3Ql l33] . Both are 
often combined in one objective function. Eor 
a set of correspondences AAAP measures the 
deviation from a global afhne transform and 
hence favors smooth deformations. ARAP lo¬ 
cally measures the deviation from a rigid trans¬ 
form and is related to elasticity since it mimics 
mechanical stiffness. These methods work on 
tetrahedral meshes, i.e., they need volumetric 
meshes that represent the shape. In the con¬ 
text of surfaces recent ARAP-like methods have 
been introduced in [ 171132 ]. 

Elastically deformable models in graphics 
have been pioneered in the end of the 1980s 
in m- Generally speaking, these models use 
characterizing properties of geometric objects 
(curves, surface, volumes) to define functionals 
that penalize deviation from them. A curve in 
three dimensions, for example, is fully charac¬ 
terized (up to rigid motion) by two parameters, 
its curvature and its torsion and hence a motion 
of the curve that changes one of them is con¬ 
sidered a deformation of the curve. Eor a vol¬ 
ume it is the local isotropic change of lengths 
that characterizes a deformation (up to rigid 
motion). This principle is reflected in the physi- 
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cal theory of (nonlinear) elasticity [HIITH]. In the 
context of shape matching deformation models 
that borrow ideas from elasticity can be found 
in [23], where the authors use linear elasticity 
for the registration of three-dimensional medical 
images, as well as in m for fitting solid meshes 
to animated surfaces. In m a thin-plate spline 
regularizer is used for point matching. Non¬ 
linear elasticity is employed in m for medical 
imaging and in [43| where two-dimensional elas¬ 
ticity is employed for three-dimensional surface 
matching via mesh parametrization. However, 
there is a vast amount of (physical) deforma¬ 
tion methods and we do not aim at discussing 
them in detail. An overview of methods used 
in graphics can be found in [38| and techniques 
used in vision and medical imaging in m- 

The second ingredient for a shape matching 
method is a data term that drives the defor¬ 
mation and/or describes dissimilarity between 
source and target shape. One way is to ran¬ 
domly sample candidate correspondences and 
to choose a deformation that best aligns the 
data. These RANSAC methods, which were 
first introduced in m, are practical for low¬ 
dimensional mapping spaces. Recently they 
were extended to deal with deformation spaces 
of higher dimension m 

Another popular group of methods comprises 
variants of the famous iterative closest point 
(ICP) algorithm, first introduced in [Hl[l5]. This 
iterative approach was first used for rigid match¬ 
ing and is well suited for various representations 
of geometric data. It is based on the computa¬ 
tion of closest points between source and target. 
Since then efficient methods have been devel¬ 
oped that use different sampling of correspon¬ 
dences, different weighting schemes that deter¬ 
mine the confidence of a candidate match etc., 
see isa [IS EO]- Recently, ICP has been ap¬ 
plied as well in the context of non-rigid match¬ 
ing O [121 HU [2H]- Some of these papers deal 
with deformations of significant magnitude. 

Techniques like multi-dimensional scaling 
(MDS) and generalized multi-dimensional scal¬ 
ing (GMDS) aim to find correspondences by 
regarding the shapes as Riemannian manifolds 
and try to find an embedding in a common met¬ 
ric space [H US]. Correspondences are then 
found in the common embedding space (which 
can even be a high-dimensional Riemannian 
manifold in the case of GMDS) rather than be¬ 


tween the shapes themselves. This way ICP can 
be regarded as an instance of the Euclidean iso¬ 
metric matching problem. MDS and GMDS are 
elegant methods but computationally costly. 

We note, however, that the vast literature and 
number of techniques can be classified also by 
other criteria. An excellent overview of the field 
and of classification criteria can be found in the 
survey paper [29]. 

1.3.2 Comparison 

This work extends our work in [46] to three di¬ 
mensions and large deformations. By favoring 
sparse forces as explanations for elastic deforma¬ 
tions we use the same prior that we used in [46] . 
There we deal with small planar deformations. 
The forces there are invariant under infinitesi¬ 
mal rotations only, due to the use of linear elas¬ 
ticity and are, in particular, not invariant under 
general rotations. 

In [16] the authors employ an enhanced ver¬ 
sion of linear elasticity by the use of rotation 
compensation, similar to co-rotated linear elas¬ 
ticity. This way they can deal with large rota¬ 
tions but this requires a large number of low di¬ 
mensional singular value decompositions in each 
iteration, whereas we solely update the tangent 
stiffness matrix computed by the FEM. Another 
difference is that in our method the optimized 
boundary forces can be made rotation invariant 
once a matching deformation is found, and it 
is theoretically suitable for very large deforma¬ 
tions. 

The splitting of the tangent stiffness matrix 
that we employ to connect boundary forces and 
boundary displacements was previously used 
in El in the context of surgery simulation and 
in [46] for matching. The technique that we 
propose generalizes this idea to the nonlinear 
case and allows us to reduce the (nonlinear) 
matching problem to solving a sequence of con¬ 
vex problems that we formulate as second or¬ 
der cone problems while simultaneously reduc¬ 
ing the computational complexity. 

In our algorithm, the computational complex¬ 
ity, i.e., the degrees of freedom to be optimized, 
is furthermore reduced by measuring the defor¬ 
mation cost, i.e., the elastic forces, on a (coarse) 
tetrahedral mesh that is needed by the FEM. 
The corresponding deformation is then extended 
to the (fine) triangular surface mesh. Similar 
ideas have used in [3Ql [33i [48] . Furthermore, 
changing the coarseness of the tetrahedral mesh 
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gives the user control over the computational 
complexity of the FEM part of the model while 
the coarseness of the surface mesh is fixed. 

The matching and the computation of the ini¬ 
tial guess is done on the fine scale of the sur¬ 
faces. We project the source points onto the 
target favoring similar points (descriptor aided) 
during the first few iterations and then switch 
to Euclidean projection once a better alignment 
is reached. This way one can regard our method 
as a two-scale version of a non-rigid ICP-like al¬ 
gorithm since the matching is done on the fine 
scale while the cost, i.e., the boundary forces, is 
measured on the coarse scale of the volumetric 
mesh. 

The paper is organized as follows. In Sectionj^ 
we give a concise overview of the model of a hy¬ 
perelastic body and describe a nonlinear EEM 
that is used in the optimization procedure. The 
optimization is described in detail in Section 
Experiments, implementation details and a de¬ 
scription of the two-scale algorithm are given in 
SectionlU Section [5] concludes with a discussion. 

2 Model and Numerical 
Methods 

2.1 Hyperelasticity 

The purpose of this section is to introduce the 
reader into the most important principles of 
nonlinear elasticity theory that are relevant in 
the sequel. As mentioned, our model optimizes 
elastic forces and hence it is important to un¬ 
derstand what forces we refer to and how they 
act. Nonlinear elasticity is a well developed and 
complex field. Mathematically oriented intro¬ 
ductions can be found in laiiHiiiu] and an en¬ 
gineering oriented introduction in m- 

Roughly speaking an elastic body is an open 
connected subset S C that reacts to applied 
forces with a deformation. It returns to its rest 
state (also called reference configuration) after 
the forces are removed and does not memorize 
previous deformations. The main equations of 
nonlinear elasticity describe a balance between 
the applied forces that cause the deformation 
and the internal force distribution inside the de¬ 
formed body (stress) in the spirit of Newton’s 
second law. We will elaborate on this in the 
following. 


Stress, Forces and Equilibrium. We as¬ 
sume that we are given an elastic body in its 
rest pose S C i.e., the rest pose is the vol¬ 
ume that is occupied by the body when no forces 
are applied to it. A deformation of the body is 
described by a locally injective and smooth vec¬ 
tor field ^ : S ^ i.e., any point x e B in 

the reference configuration has a corresponding 
point x^ := ^(x) e B^ := ^(S) in the deformed 
configuration. Any deformation should satisfy 
det(V^) > 0, i.e., subsets of B of positive mea¬ 
sure are mapped to subsets of B^ of positive 
measure. 

Suppose now that the rest state B is subject 
to forces that describe the action of the out¬ 
side world on it. These can either be volumetric 
forces f^, i.e., they are measured by per unit 
volume of the deformed configuration, or they 
can be surface forces measured per unit sur¬ 
face area acting on dB^. Examples for volu¬ 
metric forces are gravity or electric field forces, 
whereas examples for surface forces are pres¬ 
sure or spring forces. The resulting deforma¬ 
tion ^ induces an internal force distribution 
: B^ X E>‘^ ^ inside B^ that balances 
the external forces, the so-called Cauchy trac¬ 
tion vector. This vector field depends on two 
arguments since it measures the force per unit 
area acting through any cross section of B^ de¬ 
scribed by its unit normal G at x^ e B^. 
This can be expressed through a linear relation 
t{x^,n^) = T^{x^)n^ where : B^ ^ 
is the Cauchy stress tensor. Its rows represent 
the three tractions (normal and shear stresses) 
in three coordinate planes which are usually the 
three orthogonal canonical planes, see Eigure[2 
for an illustration. In particular, if x^ G dB^ 
and G is its corresponding outward unit 
normal, then t^(x^,n^) is the force acting 
on the boundary of B^ at this specific point. 
Boundary forces will be of importance in the 
sequel for the formulation of our optimization 
problem. 

The basic equations governing elasticity the¬ 
ory can now be formulated as follows (Cauchy 
principle): Let the volumetric forces be denoted 
by : B^ and let V C be an arbi¬ 

trary subvolume. Then 



This equilibrium of external forces and internal 
stress is an expression of Newton’s second law 
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Figure 1: Left: an illustration of a surface force at some point in an arbitrary cross section with normal 

of a body . Right: illustration of the three components of the Cauchy stress tensor a = with respect to the 
canonical coordinate planes. The normal stress is orthogonal to the considered plane shear stresses lie within the plane. 


and holds in the deformed configuration . A 
similar principle for the angular moments shows 
in addition that must be symmetric. Using 
the above and the divergence theorem we get 


div*r‘*’(:c*)+ /■*’(:£;■*’)) dre* =0. (2) 

Since this holds for all subvolumes V the equi¬ 
librium equations can be written in differential 
form as 



- div* T*(rr*) = /*(rr*) , rr* € S* , 
r*(rr*)n* = ( 7 *(rc*) , re* G 


( 3 ) 


where 5 * is a boundary force as described above. 


Unfortunately, this is not practical yet since 
equation (§ holds in the deformed configuration 
which is unknown. This can be remedied by a 
pullback to the reference configuration B. The 
main tool for this is the so called Piola transform 
given by: 

T{x) = det(V$(rc))T*(a;*)V$(rc)“^ . (4) 

It describes a transition from the Eulerian vari¬ 
able to the Lagrangian (reference) 

variable x. The quantity T is called the first 
Piola-Kirchhoff stress tensor. It measures stress 
on the deformed configuration per unit area of 
the undeformed configuration and it is not sym¬ 
metric. The Piola transform also transforms op¬ 
erators and forces in @ so that the equilibrium 
equations in the reference configuration B take 


the form 

-divT(a;) =/(a;), a; G B , 

T{x)n = g{x) , x G dB . 

The exact relations between the forces is given 
by 

/(a:) = (detV$(a;))/*(a;*), 

g{x) = (det V^(t)) \\\/^{x)~^n{x)\\ g^{x^) . 

( 6 ) 


This means that the transformed forces in gen¬ 
eral depend on the deformation <I> and its gra¬ 
dient even though this might not be the case 
for and g^. The interested reader is referred 
to [18]. The second Piola-Kirchhoff stress ten¬ 
sor, given by 

E(t)=V<I>(t)-^T(t), (7) 

is symmetric and measures stress on the unde¬ 
formed configuration per unit area of the un¬ 
deformed configuration. Equation (§ can be 
rewritten in terms of E as 

-div(V^(T)E(T)) =/(t) , xeB, 

T^{x)n = g{x) , x e dB 

where g{x) = \/^{x)~^g{x). However, we will 
use (§ as a basis for discretization and opti¬ 
mization. 


Constitutive Laws. The above given equa¬ 
tions describe equilibria of external and inter¬ 
nal forces but they do not take into account the 
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different reactions of materials to forces. Obvi¬ 
ously, steel reacts with a different deformation 
than rubber when exposed to the same external 
forces. Or, conversely, rubber and steel under¬ 
going the same deformation is caused by differ¬ 
ent forces. These material specific relations be¬ 
tween force and induced deformation are mod¬ 
eled by so-called constitutive laws and shall be 
described briefly in the following. 

Roughly speaking, a material is called elas¬ 
tic if the Cauchy stress can be written as 
a function of x G S and V^. With respect 
to (© and 0 we can call a material elastic if 
its first and second Piola-Kirchhoff stresses are 
functions of x e B and only. Such a func¬ 
tion is called a response function. This defini¬ 
tion allows a large class of response functions 
and can be restricted further by making reason¬ 
able physical assumptions. A common physical 
principle is the principle of frame indifference 
meaning that any observed quantity is indepen¬ 
dent of the observer. Furthermore, one can as¬ 
sume the independence of the material response 
on X G B. This is called homogeneity. Another 
restriction of the class of admissible response 
functions is the isotropy of the material which 
essentially means that the materials’ response 
at any point x G B is the same and is inde¬ 
pendent of the direction in which the force is 
applied. Steel is an isotropic material in con¬ 
trast to wood which deforms differently in the 
direction of its fibers than orthogonal to them. 
Combining all these assumptions one can show 
that the Cauchy stress and the second Piola- 
Kirchhoff stress are of the form 

= f{V^{x)V^{xf) (9) 

and 

I](a;) = S(Vi>(a;)^Vi>(a;)) . (10) 

Note that both functions, T and S are invariant 
under rigid motions. The first Piola-Kirchhoff 
stress, in contrast, can not be expressed in a 
form depending only on Nonetheless, 

T appears in the governing equations 0 and 
so the forces f,g, in contrast to g appearing 
in are not invariant under rotations. How¬ 
ever, with regard to 0 the difference is just a 
factor of From now on we will not distin¬ 

guish between the stress tensors T, S and their 
response functions T, S and we introduce the 
notations F = and C = F for the sake 
of brevity. Note that C, called Green strain, is 


the quantity that describes the local change of 
distances of nearby points (a Taylor expansion 
of II^11^ can easily show this). Furthermore, one 
can show that two deformations of a body that 
have the same Green strain differ from one an¬ 
other by a rigid motion only. Hence, two defor¬ 
mations of the same body can be identified if 
they have the same Green strain. 

Suppose now, that the first Piola-Kirchhoff 
stress is the derivative of a stored energy func¬ 
tion W{F), i.e., 

dW 

T{F) = —{F). (11) 

A material with this property is called hypere¬ 
lastic. Hyperelastic materials are advantageous 
for two main reasons. First, if also the forces 
in 0 are conservative, i.e., they can be written 
as the Gateaux derivative of some potential, the 
equations of equilibrium ^ can be written as a 
minimization problem. Since the forces will be 
unknown in our optimization problem this will 
be less important to us. The second advantage 
is that they allow a much more intuitive mod¬ 
eling of energetic penalties to certain kinds of 
deformations (e.g., non-isochoric deformations) 
rather than a direct modeling via the response 
function. 

The materials that we will use in the sequel 
are homogeneous and isotropic but before for¬ 
mulating their stored energy functions recall the 
principal invariants of a matrix A G 

ii(A) = tr A , 

i 2 (A) = tr Cof A , and (12) 

is (A) = det A , 

where tr A is the trace and Gof A = (det A)A~^ 
is the cofactor matrix (for invertible A). In 
terms of eigenvalues the principal invariants can 
be written as 

A (A) = Ai + A 2 + As , 

i 2 (A) = A 1 A 2 + A 2 AS + AiAs , and (13) 

^ 3 (A) = A1A2AS . 

The so-called reduced principal invariants are 
given by 

hiA) = i3iAy^\{A), 

hiA) = i 3 {A)~^^^i 2 {A) , and (14) 

JiA) = is{Ay^ . 
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For isotropic hyperelastic materials W can be 
formulated in terms of the (reduced) princi¬ 
pal invariants of the Green strain C only. We 
will make use of two hyperelastic materials in 
this work. The stored energy function of Saint 
Venant-Kirchhoff materials (SVK) is given by 

l^^svK = p{ii{C) - 3) -]-^(n(G) - 3)^ 

-|(*2(C)-3) 

= NTEf+,,tT{E^) 

(15) 


2.2 Discretization Using FEM 


Equation (19) is essentially a system of nonlin¬ 
ear partial differential equations (PDEs) with 
pure Neumann boundary conditions. Usually, 
this system is additionally supplemented with 
Dirichlet boundary conditions that prescribe the 
deformation of B on (parts of) its boundary. Eor 
now, we will not take into account the Dirich¬ 
let boundary conditions but later we will make 
a connection between the boundary forces and 
the deformation on the boundary. 


where 2E = C — I, C = F^F, and A, /i > 0 
are the Lame constants. SVK materials are 
linear material models since the second Piola- 
Kirchhoff stress is a linear function of E {C re¬ 
spectively): 

E{E) =\{iTE)lF2^iE . (16) 

A linearization of E yields the well known Hooke 
law of linear elasticity which is also called “small 
deformation-small strain setting”. Their poten¬ 
tial function is given by 

^Linear = , (17) 


System (19) is commonly discretized by 
means of the finite element method (EEM) [TOl 
[201 [3l]. Applying the EEM to nonlinear elas¬ 
ticity problems is in general not straightfor¬ 
ward. Eor incompressible or nearly incompress¬ 
ible problems, i.e., the volume change during the 
deformation is zero or very small, mixed EEMs 
need to be employed. However, in this work 
we will avoid this difficulty and allow volume 
changes. Eor the purpose of comparing shapes 
differing by an unknown deformation (due to 
unknown forces) this is a reasonable assump¬ 
tion and therefore it is possible to simply employ 
piece-wise linear EEMs. We shall elaborate on 
this in the following. 


where 2s = {E F^) — 21. SVK materials, 
however, can be shown to be first order approx¬ 
imations of isotropic materials in terms of E and 
therefore this model is called “large deformation- 
small strain”. 

Neo-Hookean materials have a stored energy 
function of the form 


EEMs rely on the variational form of the 
PDE, i.e., we multiply the first equation in (19) 
with a test function T in a space V of suit¬ 
ably smooth test functions and integrate the 
left-hand side by parts. Then one looks for a 
solution ^ G U such that 


W^NEO = a{h {C) - 3) + /3(J(C) - If (18) 

where a, ^ > 0 are positive constants and 
are used to model rubber and elastomers |2Q| . 
Sometimes different variants of the term involv¬ 
ing J are described in the literature. The use 
of the reduced invariants ensures the positivity 
of the stored energy. This model is the simplest 
model that is appropriate for large strain and 
large deformations. Using equation © can 
be rewritten for a general hyperelastic material 
as 


dW 

-div—(V$)=/(a;), 


dW 

W 


{V^)n = g{x) , 


X €B , 

X e dB. 


(19) 




VT dx = / /-Tdx 

Jb 


+ [ g^^dS VT G U. (20) 
JdB 

where A : B := tr(A^H) denotes the Erobenius 
scalar product. V is usually some Sobolev space 
that encodes the Dirichlet boundary conditions 
and regularity requirements on V^. We will, 
however, not get into the details of the analytical 


treatment of (20) which can be found in [MISI]. 


In order to find an approximate solution 
EEMs replace the function space U by a finite 
dimensional subspace V^. The function space 
includes an approximation B^ of the domain 
B. The problem is then to find G such 
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that 


tions of the form 


/ / /•«''* da; 

Jjsh ob J^h 

+ [ g-^'^dS G y'*. (21) 

JdB^ 

This problem is nonlinear in the deformation 
and can be solved given appropriate boundary 
conditions by Newton-Raphson type methods 
which require linearization. Let ^ t)e a 
given deformation. Then, the linearized version 
of (21) around given by 






[ f-^’^dx+l . 

JB^ JdB^ 




( 22 ) 


Note that the second derivative of W with re¬ 
spect to the deformation gradient is a tensor of 
fourth order. A system of this type has to be 
solved in each step of a Newton-Raphson proce¬ 
dure. 


Next, in order to get to a more concrete repre¬ 
sentation, we choose a global basis of the finite 
dimensional space at? ^^lat test¬ 

ing (22) with all functions of is equivalent to 
testing against all base functions. We make the 
ansatz 




i = = 1,2,3, (23) 


with real coefficients whereas we use Ein¬ 
stein’s sum convention. Note that the base func¬ 
tions are vector valued and I stands for the index 
of their relevant component. Let furthermore 
= ^0 be the expansion of ^^le basis 

of . Plugging all this into (22) we get 



This is a linear equation for the nodal deforma- 


= /Ad + gi^T) 

Vj = l,...,Ar, m=l,2,3. (25) 

Here corresponds to the term of zero order 
of the linearization of (21) and a$Q is a bilinear 
form, corresponding to the first order term of 
the linearization, whose representation matrix 
in terms of the FEM base functions is called the 
tangent stiffness matrix. The right-hand side 
represents a force term. 

As mentioned above, we use piece-wise linear 
finite elements to approximate For this pur¬ 
pose the approximate geometry is a tetrahe- 
dralization of B. A linear function in each tetra¬ 
hedron (tet) is then uniquely determined by its 
nodal values. A set of global base functions is 
defined by the relation 


y{xk) = Sjkc^ 


j,k = l,...,N, (26) 


where is the l-th canonical base vector, I = 
1,..., 3, and Xk are the nodes of B^. Using this 
we can characterize as 




{T G C^{B^,R^) I T is linear 

in each tet of . (27) 


This way each coefficient in equation (23) 


can be interpreted as the l-th component of the 
deformation vector at the node Xi of the tessella¬ 
tion B^. Note, that the discretization introduces 
two errors, the error due to the approximation 
of ^ G V, and an error that is caused by the 
approximation of the geometry B. 


3 Optimization of Boundary 
Conditions 

In this section we will give a detailed descrip¬ 
tion of the optimization procedure that we use 
to tackle the matching problem. We will show 
how to approximate this nonlinear problem by 
a sequence of convex problems that can be put 
in the form of a second order cone program 
(SOCP). 

Our goal is to match two two-dimensional sur¬ 
face meshes, the source S and the target T, de¬ 
viating from one another by an unknown defor¬ 
mation Very often these surfaces represent 
the surface of actual physical entities like solid 
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bodies. In order to deform solid bodies an ex¬ 
ternal force is necessary. As stated in the intro¬ 
duction, the external forces we are seeking are 
simple, i.e., sparse and isotropic and act on the 
boundary of the physical body B only. These 
forces are a priori unknown as well as the de¬ 
formation of the surface itself. But we can put 
them in relation using the FEM, as we will show 
in the sequel. The optimization problem of find¬ 
ing the deformation can be formulated as follows 


minimize 


subject to 


r 

dW, , 

JdB 



dS 


dW 

div—(V$)=0, x€B 

Or 


^(dB) = r, 


we get a deformation of S that is caused by 
sparse forces and is “closer” to the target and 
we can re-iterate the process. In each step we 
have to solve a nonlinear optimization problem 
which can be formulated as 



/ i 

dW 

minimize 

JdB \ 

^(V$)n 


f dW 

subject to / ——(V<F) : VTdx = 

Jb 

f dW 

/ ^(V^)n.Td5 

JdB oh 

VT G E. 


(29) 


i.e., we are imposing elastic properties on the 
source S = dB and then we are seeking an elas¬ 
tic deformation ^ of S such that the external 
forces that are responsible for the deformation 
of S into T are minimal. The objective in (28) is 
essentially an L^-cost on the force. Note that we 
assume a pure boundary force and no volumet¬ 
ric forces. The force term we optimize in the ob¬ 
jective function of (28) is measured by the first 
Piola-Kirchhoff stress tensor, see equation 
They are acting on the deformed configuration 
but are measured per unit area in the reference 
configuration. Hence, the objective is not in¬ 
variant under rotation. This, however, is not a 
limitation because one can simply take the re¬ 
sult of the optimization (28) and apply transfor¬ 
mation 0 to get rotation invariant forces that 
are measured and acting on the reference con¬ 
figuration. This way one can identify two de¬ 
formations that deviate from one another by a 
Euclidean transform. 


Unfortunately, the optimization problem ( [2^ 
is nonlinear due to the second constraint and 
due to the nonlinearity of material models. We 
will show a way to approximate it by a sequence 
of simpler and, most importantly, convex prob¬ 
lems. 


In contrast to the boundary forces one can 
often give a rough guess <I>* of surface corre¬ 
spondences between the surface meshes to be 
compared. Again, regarding the source mesh S 
as the boundary of a physical body this guess 
of correspondences induces a boundary force on 
dB. This boundary correspondence can then be 
optimized according to our prior that the in¬ 
duced forces are sparse and isotropic. This way 


This way we eliminate the second constraint 


in (28) by localizing the optimization around 


the initial guess The new term in (28) can 


be interpreted as a spring force with spring con¬ 
stant k, measuring deviation from the estimated 
correspondences, that needs to be balanced by 
sparse elastic forces. Note that we also replaced 


the PDE constraint in (28) by its variational 


form (20). 


Problem (29) can be given a discrete formu¬ 
lation using equation ( [21] ), i.e.. 


minimize 


r / 

dW, 

IdB^^ Y 



+ ^ 11$'* ll^ 'l d5 


f dW 

subject to / — (V$'*) : V$'* d:c = 

J]gh oF 

f dW u u 
JdBh dh 


V G W . 


(30) 


This is still a difficult nonlinear problem al¬ 
though now finite dimensional. The constraint 
as an equation itself, i.e., the discretized vari¬ 
ational form of the nonlinear elasticity equa¬ 
tions, is usually solved by means of a Newton- 
Raphson type procedure that relies on succes¬ 
sive linearization of the nonlinear terms. This 
linearization suggests a way to treat the bound¬ 


ary force term in the objective of (30). In each 


Newton-Raphson step we have to solve a linear 


system of the type (25). This system can be re- 
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cast by a renumbering of the nodal deformations 
into the form 


Gb 


Bb,^^ 

0 


Bp^h 


A 




A 




A 


II Ao 


^B 


i.e., we decompose the deformation vector 


(31) 


(32) 


into a boundary deformation part ^b collecting 
all the deformation vectors on boundary nodes 
of and into a part of deformation vectors 
at the inner nodes. The right-hand side rep¬ 
resents external boundary forces Gb and zero 


external volumetric forces in (25). 


Equation (31) can be used to replace the non¬ 


linear constraint in (30): 


minimize 

^ IdB^ 


■I 

Jdl 


subject to 


A 

A 




- 11 ^^ - ^^’*1 

2 II I 

^lAo 
A 


:ds 


Gb 


0 



IBAo 


BIAo 

hlAo 


^B 


(33) 


problem as 


K 


minimize 


E r 


i=l 


^ + (37) 


Each GB,i, i = 
attached to the 


which is a convex problem. 

1,..., K, is a force vector in 1 
i-th of the K boundary nodes. Note that (37) 
is unconstrained and that GB,i depends on^^g 
and linearly on Also, the dimension of the 
(local) optimization problem is significantly re¬ 
duced since it only involves the boundary defor¬ 
mations ^B • Such condensation techniques have 
been used in case of completely linear elasticity 
in [TTlHB] . 

The derived local optimization problem 
describes a single step in our optimization 
scheme and needs three initialization parame¬ 
ters: 


• An initial deformation needs 

to be given in order to assemble the tangent 
stiffness matrix. We set = I before the 
first step. 


• The spring constant k determining the 
strength of the spring force. As the source 
gradually deforms into the target during 
the iterations the elastic force that resists 
the springs increases. Thus, one needs to 
increase k as well. 


Note that the boundary force term Gb in the 
linearized constraint is essentially a smoothed 
version of the boundary force term 


dW 

~W 




(34) 


in the objective since it is weighted by the finite 
element base functions. 


Equation (31) is in a form such that it is pos¬ 
sible to eliminate the inner deformations by 
taking the Schur complement with respect to the 
block Abb 


Gb — BB^^h — Abi^^hA^^ 


+ S^h 


where 


(35) 


• An initial guess of correspondence : 

S ^ T determining the direction of the 
spring force. This initial guess becomes 
more accurate the “closer” the source is to 
the target shape. 


The last point suggests that it makes sense 
to constrain the search for improved correspon¬ 
dences locally to the target shape T. This is 
done by increasing the spring force penalty for 
deviation from the estimated correspondences 
into normal direction during the process and 
amounts to replacing the Euclidean metric in 


the second term of (37), i.e.. 


K 


- mil = N ii^B(^i) - m^i)\\2 (38) 


i=l 


S^h — Abb,^^ — ABi^^hAj^^^hf^jB^^h . (36) 


Utilizing (35) we can rewrite the optimization 


with a Mahalanobis distance that depends on 
the corresponding point <l>^(x^) G T where Xi is 
a node of S. Such a local metric is described by 


10 

































July 29, 2015 


Nonlinear Force Optimization for Shape Matching 


a symmetric positive definite 3-by-3 matrix 




■ ■ 




. ^2 . 



X. 


0 0 
0 A, 0 

0 0 A, 


(39) 

where n is the unit normal at G T and 

ti,t 2 form an orthonormal basis of the tangent 


space. We can then write (38) as 

K 


where 


i=l 


(40) 

(41) 


bllg = {x,Qx) 

and Q = diag((3i,..., Qk)- Decreasing the ra¬ 
tio Xn/Xt then means increasing the penalty in 
normal direction. Combining (40) with (37) and 


introducing slack variables gi and e we can re¬ 


formulate (37) as 


K 

minimize Oj 
subject to Gb.i 


k 

2 ® 


< ft ,* = 1, 


.K 


(42) 




The first K constraints are already in second 
order cone form and the constraint on e can be 
reformulated using 

\\^B-nfQ<\{ie + lf-{e-lf) (43) 


into 


2R{^b - 
e — 1 


< e 


(44) 


where R = diag((5i,. •., ^ The SOCP in 

each iteration then becomes 


. . . ^ k 

minimize > gi -\—e 




subject to 


Gs.i 


<ft = 


(45) 


2i?($B - 

e — 1 


< e . 


There are several efficient solvers for SOCPs. 
We found that MOSEK [4] is very suitable and 
used it through the YALMIP interface [37] with 
MATLAB. 


4 Implementation and 
Results 


[ n ti t 2 ] 


4.1 The Two-scale Algorithm 

In this section we describe our algorithm and 
its implementation aspects in detail. As in¬ 
put our method takes two surface meshes S, T. 
The source S is assumed to represent the sur¬ 
face of a deformable body. Elastic properties on 
S can then be imposed on a volumetric tetra- 
hedralization . Our formulation so far as¬ 
sumes that the boundary of the tetrahedraliza- 
tion coincides with the surface mesh. Tetrahe¬ 
dral meshes that respect the given surface tes¬ 
sellation S can be obtained, for example, with 
TetGen |42) . This, however, has two disadvan¬ 
tages. Eirst, the mesh S is required to be free of 
artifacts like self-intersecting triangles or non¬ 
manifold edges. This is, in general, not the case 
for meshes obtained from scanned data. Sec¬ 
ondly, a tet mesh with the same boundary as 
S will have a large number of degrees of free¬ 
dom that are relevant in our model and although 
the optimization procedure is effectively reduced 
to the boundary of this is computationally 
very expensive. We can remedy this by measur¬ 
ing elastic properties on a coarser scale than the 


spring force in (37), as we explain below. 


To this end we created a coarse tetrahedral- 
ization B^, typically with 2000 to 3000 tets 
(compared to the surface meshes that contained 
between 10000 to 50000 triangles). We used the 
Iso2Mesh toolbox [22| that wraps code provided 
in m and [42|. Next, we extracted the bound¬ 
ary mesh of B^ and projected every node in S 
onto the mesh dB^, i.e., for each node p e S 
we computed the triangle t G dB^ with minu- 
mum distance and the corresponding projection 
point Xt ^ t. We then computed its barycentric 
coordinates {u, v, w) in t. The deformation of p, 
4>(p), can then be written as a linear combina¬ 
tion 

$ (p) = uifg (a:;i) + v(fi% (x2 )+w<p^(x3) (46) 

of the deformation of the nodes xi,X 2 ,X 3 of t. 
Erom now on we will therefore distinguish be¬ 
tween the deformation <I> : 5 ^ and the 
deformation : dB^ of the boundary 


of the tet mesh. An equation of the form (46) 


applies to each displacement vector 4>(p) for any 
p G 5. Therefore, any deformation of dB^ 
can be mapped linearly to a deformation of S. 
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deformed 
and forces 


Algorithm 1 : Sketch of the algorithm 
input : Surface meshes S, T and elasticity 
parameters 

output: Deformations 
meshes ^{S), 

Gb 

1 begin 

2 Create a tetrahedralization of S; 

3 Project the nodes of S onto the mesh 

dB^ and compute the smoothing map 
T : C^{dB^,R^) described 

in (|47|); 


9 

10 

11 


13 


14 

15 

16 


17 


18 

19 

20 


Precompute a local orthonormal 
coordinate system at each point g G T 
to set up a local metric as described 
in (39); 

Set*^ = I, ^ = I, ^ = 0 and k = 1; 
while k < Kmax E is not 
stationary do 

Compute an initial guess of 
correspondences : ^{S) —> T; 
For each point chose 

the corresponding metric QJ. as 
in (39); 

SetV= 1; 

while I < Lmax do 

Compute the first order 


approximation (31), i.e., the zero 


order term and the tangent 
stiffness matrix at ip^', 

Solve the local optimization 


problem (48) with initial 
deformation in (49) to obtain 

Use the optimized boundary 
conditions PB^k • ^ 

solve for the full deformation 
p^ : B^ ^ R^ by means of a 
Newton-Raphson solver; 

Update Pq = p^ and / = / + !; 
end 

Save the stored energy E = p^Gb of 
the deformed body and optionally 
various other elasticity measures 
(e.g., strain, stress, etc.); 


Increase the spring constant in (48) 
and possibly decrease the ratio 
An/At in (1^; 


Update k ^k + 1; 


end 

Return; 


21 end 


In addition we smooth the interpolated defor¬ 
mation field with a number of damped Jacobi 
relaxations using the Laplace-Beltrami opera¬ 
tor A^. Note that smoothing the (linearly in¬ 
terpolated) deformation field on S amounts to 
smoothing the mesh 4>(5). Interpolation and re¬ 
laxation can be represented by a linear operator 
T that acts on the space of continuous functions 
G^{dB^,R^) and maps into (7^(5,M^), i.e., we 
can write 

$ = tPb ■ (47) 


We note that an alternative approach to ex¬ 
tending a deformation from a coarse to a fine 
scale was taken by Kovalsky et al. m and 
in [48]. There the map T was determined by 
a moving least squares optimization, i.e., each 
deformation vector 4>(p) on S is the weighted 
average of nearby deformation vectors on B^. 
This way local volumetric deformation of B^ in¬ 
duces a local deformation on S whereas we use a 
surface deformation to surface deformation ap¬ 
proach. 


We can now reformulate problem (45) as fol¬ 
lows. 


K 


E rv 

Qi H—^ 


minimize 

gi,e,v^,4> ^ 

subject to K?B,i 


< Qi = 


- $*) 
e — 1 


< e 




(48) 


where 

Gb + S^hp^ 

(49) 

in the spirit of ( [3^ . This is an optimization 
problem on two scales since it measures the 
(computationally expensive) deformation cost 
for the volumetric body on a coarse scale while 
the spring force between source and target shape 
is measured on a fine scale. The algorithm is 
summarized in Algorithmic 

Our Matlab implementation was run on a 
standard desktop PC with 8GB of RAM and In¬ 
tel Core i7 with 3GHz. The bottleneck routines 
such as the assembly of the tangent stiffness ma¬ 
trix were written in C. For SVK materials we 
wrote our own FEM code and for Neo-Hookean 
materials we modified the source code of Cal- 
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Figure 2: Deformation tracking experiment. The deformation sequence consists of 60 frames 7fc. In each frame the upper 
face of the beam S was twisted until a 360 degree twist was reached. The source S was then matched to each frame 7k while 
storing the result and gradually increasing the deformation. We show the result of the matching for four frames 

(90°, 180°, 270° and 360° twist) for the Neo-Hookean material model (upper row) and the Saint Venant-Kirchhoff Model 
(lower row). The optimized force vectors, measured on the coarse mesh 23^, are shown on the results. The middle row 
shows the corresponding target frames (turquoise). The surface mesh S consisted of 14400 triangles while the tetrahedral 
mesh on which the forces of the deformation are measured consisted of 1536 tets (567 surface triangles on dB^). 


culiX [21], a three-dimensional FEM software 
that is available under the GNU General Public 
License. 


Computing the Initial Guess. Comput¬ 
ing a reasonable initial guess is not straightfor¬ 
ward and is mostly done in a heuristic man¬ 
ner. In the case that S and T only deviate 
by small deformations and rotations it is often 
reasonable to find initial correspondences sim¬ 
ply by projecting points of the source mesh onto 
the target as done in a standard iterative clos¬ 
est point algorithm. In the case of large defor¬ 
mations and/or large rotations, however, near¬ 
est neighbor projection is usually a bad choice. 
The limitation of small deformations and rota¬ 
tions can partially be overcome by the use of 


shape descriptor like the heat kernel signature 
(HKS) [49] or wave kernel signature (WKS) [?]• 
These descriptors describe the diffusion of heat 
at each point or the dispersion in case of the 
WKS. Since these descriptors are based on the 
spectrum of the Laplace-Beltrami operator they 
are invariant under isometries, i.e., surface de¬ 
formations that preserve geodesic distances. For 
large isometric or nearly isometric deformations 
one can then use a nearest neighbor search in 
the descriptor space for the first few iterations 
of Algorithm instead of the Euclidean nearest 
neighbors. However, the so obtained correspon¬ 
dences can be treacherous, in particular in the 
presence of inner symmetries of the shape. 

In our work we use a descriptor guided near¬ 
est neighbour search in the Euclidean space. As 
a shape descriptor we used the HKS. In the 
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first few iterations we seek for each point on the 
source shape S k-nearest neighbours and then 
we select the one with the best matching HKS. 
We further assign a confidence value w G [0,1] 
for the correspondence and use only correspon¬ 
dences with confidence above a certain thresh¬ 
old. Similar weighting procedures have been 
used in (301 HSl Ei. The number k of com¬ 
puted nearest neighbors was usually set around 
two to four percent of the number of nodes of 
the shape. After a few iterations we decreased 
k down to one and passed on the confidence 
weighting what amounts to standard Euclidean 
nearest neighbours. 


4.2 Experiments 


Material Models. For our experiments we 
used the three material models introduced in 
Section In the following we will refer to 
the small deformation model as the linear 
model, to the Saint Venant-Kirchhoff model (15) 
as SVK and to the Neo-Hookean model (18) as 
NEO. 


The first experiment is shown in Figure 
and is supposed to demonstrate the flexibility 
of the two large deformation models (SVK and 
NEO). We took a triangular surface mesh S of 
a beam and created a sequence of deformations 
Tk (which we will call frames). In each frame we 
twisted the upper face of the beam while keep¬ 
ing the bottom face fixed until a twist of 360 
degree was reached. We then used the frames 
Tk as target frames for matching S according 
to Algorithmic We start from frame k = 1 and 
gradually increase the deformation while storing 
the resulting deformations ^k{S) until we reach 
the last frame. This is a very large deformation 
that causes high mechanical stress. Note that 
all mechanical quantities are measured with re¬ 
spect to the undeformed configuration S. 

Here, we did not perform the smoothing de¬ 
scribed in step 3 in Algorithn|C f^ie first few 
steps both materials performed well. After a 
twist of approximately 270 degree we observed 
that the SVK model needed more iterations in 
each step to converge than the Neo-Hookean 
model. Also, the Newton-Raphson solver that 
is applied after every single optimization step 
needed more iterations to converge to an actual 
solution of the elasticity equations. The Neo- 
Hookean model, in contrast, showed a substan¬ 
tial computational robustness against this large 
deformation. The linear material model, i.e.. 



Figure 3: Our algorithm can find rigid motions. The tar¬ 
get meshes T (transparent green, right) are rotations of the 
source meshes S (left). The bird mesh was taken from the 
SHREC data set. Our tet mesh had approximately 1400 tets 
compared to 18000 triangles of the surface mesh. 



Figure 4: Discrete L^-norm of the boundary forces in each 
iteration. The search space for the deformations is not con¬ 
fined to rigid motions but the deformation is “not far” from 
a rigid motion. 


the SVK model with linearized strain (no up¬ 
date of the tangent stiffness matrix in step 11 
of Algorithmic failed to produce reasonable re¬ 
sults much earlier. However, we do not show the 
linear case here since this material model is not 
suitable for this large deformation anyway. 

This coincides with our expectation since the 
SVK model is appropriate for large deforma¬ 
tions but not suitable for large strain in contrast 
to the Neo-Hookean model which is suitable for 
both. 

Invariances. As elaborated in Section |2] 
for each nonlinear material the deformation 
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Figure 5: Several steps of the matching of a source shape (left, gray) to two different articulated human shapes (right, 
green) from the SCAPE dataset are shown in the lower and upper row. The approximating volumetric mesh is shown 
on the left (orange) as well as the force vectors measured on the volumetric mesh. The force vectors are scaled for better 
visibility. The surfaces have approximately 25000 triangles. Our volumetric tessellation had approximately 2000 tets. 


force measured by the second Piola-Kirchhoff 
stress 0 is invariant under rotations. How¬ 
ever, in our numerical scheme we optimize a 
smoothed version of forces derived from the 
first Piola-Kirchhoff stress which is, in general, 
not invariant under rotations. Nevertheless, the 
magnitude of the forces is invariant but not the 
direction. For pure rotations the force is zero. 
Hence, we would like the matching algorithm to 
find pure rotations and to return a zero force. 
Figure shows two examples. The initial guess 
of correspondences for the bird was made and 
updated using a nearest neighbour search in de¬ 
scriptor space in the first few iterations of our 
algorithm. For the rotated beam we simply used 
Euclidean projections since the HKS is “less de¬ 
scriptive” due to the high number of symmetries 
and hence might lead to an unreasonable initial 
guess. 

The graphs in Figure]^ show the norm of the 
optimized boundary force in each iteration of 
the optimization. One can clearly see that in 
the first few iterations the optimized boundary 
force is not zero. This is to be expected since 
we do not restrict the search space of deforma¬ 


tions to rotations. In the last few iterations the 
force is close to zero implying that our algorithm 
found a rigid motion of The slightly incor¬ 
rect matching of the rotated surface mesh of the 
bird, as seen in Figure]^ (bottom right), is due 
to the interpolation of the deformation of B^ 
to S and can be remedied by taking a finer tet 
mesh at the cost of computation time. 

Last, we should mention that in case of a 
linear material model (Hooke’s law) the forces 
are not invariant under rotations but under in¬ 
finitesimal rotations. These maps are the gen¬ 
erators of the rotation group, i.e., they form the 
Lie algebra of 5'0(3), and have skew symmetric 
derivatives. In two dimensions the effects of this 
invariance were demonstrated in [46]. Here, we 
will refrain from showing this. 

More Experiments. Figureshows two ex¬ 
amples of matching human poses. The models 
were taken from the SCAPE dataset [5]. Both 
poses were aligned using the SVK model. Note 
that the SCAPE dataset is nearly isometric and 
hence the HKS is nearly invariant under change 
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Figure 6: Examples from the SHREC dataset 1261 . The source shapes are shown in the top row. We compare three 
material models: the linear model, the Saint Venant-Kirchhoff model (SVK) and the Neo-Hookean model (Neo). Several 
manually labeled landmark points of the source meshes are shown on the deformed source mesh and then compared to 
their corresponding landmarks in the target shapes (green). The SHREC surface meshes we used consist of 9000 to 14000 
triangles and the tetrahedral meshes of 1600 to 2700 tets. 


of poses. This facilitates the nearest neighbor 
search in the first few iterations of our algo¬ 
rithm. Note that the external force found by 
our algorithm is sparse and mostly acts on the 
articulated parts of the model. 

Examples of matching taken from the 
SCHREC dataset [26j are shown in Eigure 
Note that this dataset consists of usually non¬ 
isometric instances within each object class. 
Still, our algorithm produces good results. 

Two examples taken from the TOSCA 
dataset [13] (left and right poses) are shown 
in Eigure [7 from two different viewpoints (top 
and bottom). We used the SVK and the linear 
model, i.e., the SVK model without update of 
the tangent stiffness matrix (Hooke’s law), for 
each pose. Both results exhibit similar quality 
whereas the SVK model exhibits deformation 
forces that are higher in magnitude. This is in 
line with the comparison of the SVK model and 
the linear model in the experiment shown in Eig¬ 
ure [6l 

The experiments performed so far have been 
carried out on closed surfaces without boundary. 



Linear 

SVK 

NEO 

Birds 




Norm of forces 

6 

14 

2 

^ Iterations 

14 

35 

36 

Planes 




Norm of forces 

6.01 

9.4 

7.9 

^ Iterations 

17 

21 

36 

Pliers 




Norm of forces 

0.6 

4.7 

0.8 

# Iterations 

7 

7 

11 


Table 1: Comparison of the three material models used on 
the examples of Figure We show the sum of the norms of 
the force vectors and the number of iterations 


However, the method is not restricted to closed 
surfaces provided one can create a volumetric 
mesh that represents the given object. To create 
the volumetric mesh it is usually necessary to fill 
holes. This, however, is the only step at which 
we use the closed surface. The rest of the steps 
follows the steps described in Algorithm We 
demonstrate this in Eigure There we took 
two models of teeth from the anatomical dataset 
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Figure 7: Two examples from the TOSCA dataset 1131 from two different viewpoints (upper and lower row). The source 
(middle) is matched to the targets (green). Each shape in the TOSCA dataset consists of approximately 55000 triangles. 
The volumetric mesh we used consists of 1500 tets. Both results are of similar quality whereas the SVK model needed 
forces of higher magnitude. 


of [9] which exhibits quite large deformations. 
The deformation of the underlying tetrahedral 
mesh is shown in Figure 

5 Discussion 

In this work we propose a new ICP-like method 
in combination with a nonlinear regularizer for 
three-dimensional surface matching. The idea is 
to regard the surface as the surface of a physical 
body that can undergo elastic deformation. Our 
approach to finding a reasonable deformation is 
led by the assumption that the forces acting on 
the body are of a simple nature, i.e., we assume 
the forces are sparse, isotropic and act on the 
boundary only. This is, we believe, a reasonable 
assumption in many scenarios of shape deforma¬ 
tion such as the change of a human pose. 

We model the physical body as a hyperelas¬ 
tic material. The underlying equations, a non- 
linearly coupled system of PDFs, is solved by 
means of a conformal FEM and needs to be 
supplemented with appropriate boundary condi¬ 
tions. We provide pure Dirichlet boundary con¬ 


ditions, i.e., the boundary deformation. This 
boundary deformation is in general unknown 
but one can give an initial guess. This initial 
guess induces a spring force between the source 
and the target that is supposed to drive the de¬ 
formation of the source. We then use our pro¬ 
posed method to improve the initial guess such 
that the forces that act on the deformed body 
are sparse and isotropic according to the prior. 
This is done using the FEM which allows us to 
split the tangent stiffness matrix. The splitting 
we employ allows to connect boundary forces on 
the deformed body and boundary deformation 
of the source. This splitting generalizes the lin¬ 
ear condensation methods that have been used 
in [nmi to the nonlinear case. 

The optimization problem to be solved is non¬ 
linear and computationally complex, depending 
on the shape and on the material model. We 
show how to approximate this problem by a 
sequence of convex problems that is solved in 
an iterative fashion with low computational de¬ 
mands. The convexification is done by succes¬ 
sive linearization of the material model. The 
computational complexity is lowered by approx- 
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Figure 8: Two experiments performed on tooth shapes taken from the dataset of anatomical shapes of (top and bottom 
row). The meshes have disk topology and the deformations are nonisometric and exhibit large strain. Several intermediate 
steps of the matching procedure of the sources (left, gray) into the target (right, green) are shown. We used the Neo- 
Hookean model. Each mesh we used has about 10000 triangles. The volumentric meshes that we created had about 1000 
tets. 


imating the source surface with a coarser tetra¬ 
hedral mesh. A deformation of the boundary 
of the volumetric tetrahedral mesh then natu¬ 
rally induces a deformation of the source surface 
by means of interpolation and smoothing. This 
amounts to a two-scale algorithm since the elas¬ 
tic properties are measured on the coarse scale 
while the provided initial guess is a map on the 



Figure 9: The volumetric mesh in rest state approximating 
the source shape shown in the top row of Figure is shown 
on the left. The right side shows the deformed mesh after 
matching including the conformal distortion of each tet. 


(fine) source surface. Reducing the optimization 
problem to the boundary of the coarse scale vol¬ 
umetric mesh further reduces the dimensionality 
of the problem significantly. 

Our method often delivers good results but 


can fail in certain cases. The estimated corre¬ 
spondences that need to be given before each it¬ 
eration of the optimization is, at least in the first 
few iterations of our algorithm, usually com¬ 
puted by a Euclidean /^-nearest neighbor search. 
Out of these k nearest neighbors we pick the one 
with best matching HKS. In case of strongly 
non-isometric deformations or shapes with a 
high amount of intrinsic symmetry this will re¬ 
sult in very unreasonable initial correspondences 
and the algorithm will get stuck in a local min¬ 
imum. The case k = 1 is the usual Euclidean 
nearest neighbor search and will give unreason¬ 
able correspondences in case of large deforma¬ 
tions and/or rotations even if the deformation 
is isometric. 

Another drawback is due to the coarse vol¬ 
umetric meshes. Eor example, a pure rotation 
of a coarse mesh does not necessarily induce a 
pure rotation of a fine mesh. The reason for 
this is that we essentially map a linearly inter¬ 
polated map of the coarse mesh to the fine scale 
as described in Section [4T] without incorporat¬ 
ing information in normal direction. This effect 
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is shown in Figure Also, smoothing the de¬ 
formation after interpolation can smear features 
such as corners of the deformed surface in an 
undesired way. Taking a finer volumetric mesh 
reduces this effect. 



Figure 10: The effect of the interpolation of a pure rotation 
of the coarse mesh on the fine surface mesh. We chose a 
rotation R E_SO(3) and constructed the map T according 
to equation The rotated coarse mesh (orange) is shown 

in the upper row. The bottom row shows the effect of the 
interpolated rotation on the fine mesh. We compare this to 
the rotated fine mesh (turquoise). 

Furthermore, the method is not guaranteed to 
avoid element flipping due to the linearization 
and can exhibit high conformal element distor¬ 
tion, see Figure Methods dealing with this 
problem in two and three dimensions can be 
found in 1301 El im and can be integrated into 
our framework. 

Finally, the computed boundary forces that 
are responsible for the shape change can be 
made invariant under rigid motion by a pull¬ 
back to the undeformed configuration. Hence, 
we believe, that they can serve as a comparison 
measure for different deformations of the same 
object since they measure how “difficult” it is to 
achieve the observed change. 
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